home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
a_utils
/
expanded.lha
/
expanded
/
src
/
Basis.C
next >
Wrap
C/C++ Source or Header
|
1992-03-19
|
27KB
|
1,029 lines
//
// Linear-Affine-Projective Geometry Package
//
// Basis.C
//
// $Header$
//
// William J.R. Longabaugh
// University of Washington
//
// Implementation of the linear-affine-projective geometry
// package described in William J.R. Longabaugh, "An Expanded
// System for Coordinate-Free Geometric Programming", Master's
// thesis, University of Washington, 1992.
//
// Copyright (c) 1992, William J.R. Longabaugh
// Copying, use and development for non-commercial purposes permitted.
// All rights for commercial use reserved.
// This software is unsupported and without warranty; it is
// provided "as is".
//
// ***********************************************************************
#include <math.h>
#include "Lap.h"
// ***********************************************************************
// ***********************************************************************
//
// Basis Class
//
// ***********************************************************************
// ***********************************************************************
//
// Private/protected member functions
//
// ***********************************************************************
//
// Used by other bases to build new bases:
//
Basis::Basis(BasisType b, Space& ins,
char *namein, Matrix &m) : s(ins), tostdb(m), fromstdb(Inverse(m))
{
type = ANY_BASIS;
holds = b;
strcpy(name, namein);
}
// ***********************************************************************
//
// Used to build a standard basis:
//
Basis::Basis(BasisType b,
Space& ins, int size) : s(ins),
tostdb(IdentityMatrix(size)),
fromstdb(IdentityMatrix(size))
{
type = ANY_BASIS;
holds = b;
strcpy(name, "Standard Basis");
}
// ***********************************************************************
//
// Dumps out Basis data:
//
void Basis::data_out(ostream &c, int indent, char* n)
{
char *ibloc = new char[indent + 1];
for (int i = 0; i < indent; i++) {
*(ibloc + i) = ' ';
}
*(ibloc + indent) = '\0';
c << ibloc << ast;
c << ibloc << n;
c << ibloc << "Type of basis: ";
BasisTypeOut(c, type);
c << "\n";
c << ibloc << "Currently holds: ";
BasisTypeOut(c, holds);
c << "\n";
if (holds != NULL_BASIS) {
c << ibloc << "Name of basis: " << name << "\n";
c << ibloc << "Contained in space:\n";
s.debug_out(c, indent + ERR_IND);
c << ibloc << "Matrix representation of basis wrt standard basis:\n";
tostdb.debug_out(c, indent + ERR_IND);
c << ibloc << "Matrix representation of standard basis wrt this basis:\n";
fromstdb.debug_out(c, indent + ERR_IND);
}
c << ibloc << ast;
delete ibloc;
return;
}
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Need to do memberwise initialization, since Matrix class members need
// to do memberwise initialization.
//
Basis::Basis(Basis& b) : s(b.s), tostdb(b.tostdb), fromstdb(b.fromstdb)
{
type = ANY_BASIS;
holds = b.holds;
b.Name(name);
}
// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment.
//
Basis& Basis::operator=(Basis &b)
{
//
// This unsavory checking is required to bypass the apparent inheritance of
// assignment under g++ 1.37:
//
if ((type != ANY_BASIS) &&
((type != b.holds) && (b.holds != NULL_BASIS))) {
errh.ErrorExit("Basis& Basis::operator=(Basis&)",
"Illegal assignment attempted", *this, b);
}
holds = b.holds;
s = b.s;
b.Name(name);
tostdb = b.tostdb;
fromstdb = b.fromstdb;
return (*this);
}
// ***********************************************************************
//
// Output for debugging:
//
void Basis::debug_out(ostream &c, int indent)
{
static char* name = "Basis object\n";
this->data_out(c, indent, name);
return;
}
// ***********************************************************************
//
// Return the nth member of the basis.
//
GeOb Basis::operator[](int n)
{
static char* sig = "GeOb Basis::operator[](int)";
int maxn = s.MatrixSize();
int i;
GeObType rettype;
Space retspace(s);
Matrix hold;
// First check that the integer is in the correct range. Note that
// for projective spaces, the user can request the unit point of
// the projective frame.
if (holds == HFRAME) {
maxn++;
}
if ((n >= maxn) || (n < 0)) {
errh.ErrorExit(sig, "Specified n is invalid",
ErrVal("Value of n = ", n), *this);
}
// The members of the basis are encoded in the matrix. We need to
// extract the correct row and build a GeOb using that row. Two
// special cases are the vector components of frames and the
// unit points of projective frames:
if ((holds == FRAME) && (n != maxn - 1)) {
retspace = s.GetSpace(TANGENT);
rettype = AFF_VECTOR;
hold = (tostdb * Transpose(s.HoistTangent()))[n];
} else if ((holds == HFRAME) && (n == maxn - 1)) {
rettype = PROJ_POINT;
hold = tostdb[0];
for (i = 1; i < tostdb.Rows(); i++) {
hold += tostdb[i];
}
} else {
rettype = s.NativeType();
hold = tostdb[n];
}
GeOb retval(rettype, retspace, hold);
return (retval);
}
// ***********************************************************************
//
// Apply the basis to the geometric object to get its coordinates
// with respect to the basis.
//
ScalarList Basis::operator()(GeOb &v)
{
static char* sig = "ScalarList Basis::operator()(GeOb &)";
Boolean cannot_map = FALSE;
GeOb new_obj;
Matrix retvals;
int i;
GeObType native = s.NativeType();
// Check to see if object can be cast into the appropriate object
// in the space of the basis, and do any necessary casting. Then
// get the coordinates with a matrix multiply. Special case where
// the basis is in an affine space and the object is from the
// tangent space:
if (v.CanMapTo(native)) {
new_obj = v.MapTo(native);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else {
retvals = new_obj.MatrixRep() * fromstdb;
}
} else if (((holds == SIMPLEX) || (holds == FRAME)) &&
v.CanMapTo(AFF_VECTOR)) {
new_obj = v.MapTo(AFF_VECTOR);
if (new_obj.SpaceOf() != s.GetSpace(TANGENT)) {
cannot_map = TRUE;
} else {
retvals = (new_obj.MatrixRep() * s.HoistTangent()) * fromstdb;
}
} else {
cannot_map = TRUE;
}
// Error exit if we cannot match spaces:
if (cannot_map) {
errh.ErrorExit(sig, "Argument cannot be mapped to space of basis",
v, *this);
}
// All that remains is to build a scalar list from the coordinates:
ScalarList retlist(retvals.Columns());
for (i = 0; i < retvals.Columns(); i++) {
retlist[i] = retvals[0][i];
}
return (retlist);
}
// ***********************************************************************
//
// Return the dual basis for a vector basis.
//
VBasis Basis::Dual(void)
{
// Must be a vector basis:
if (holds != VBASIS) {
errh.ErrorExit("VBasis Basis::Dual(void)",
"Basis is not a vector basis",
*this);
}
// Copy the basis name, prepending with "Dual to". If it is the
// dual of a dual, return the original name:
char buf[MAX_NAMESIZE];
char *dualto = "Dual to ";
if (strncmp(name, dualto, 8) == 0) {
strcpy(buf, name + 8);
buf[MAX_NAMESIZE - 1] = '\0';
} else {
strncpy(buf, dualto, 8);
strncpy(buf + 8, name, MAX_NAMESIZE - 8);
buf[MAX_NAMESIZE - 1] = '\0';
}
Basis retval(VBASIS, s.Dual(), buf, Transpose(Inverse(tostdb)));
return ((VBasis)retval);
}
// ***********************************************************************
// ***********************************************************************
//
// Simplex Class
//
// ***********************************************************************
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment:
//
Simplex& Simplex::operator=(Simplex &b)
{
holds = b.holds;
s = b.s;
b.Name(name);
tostdb = b.tostdb;
fromstdb = b.fromstdb;
return (*this);
}
// ***********************************************************************
//
// Output for debugging:
//
void Simplex::debug_out(ostream &c, int indent)
{
static char* name = "Simplex object\n";
this->data_out(c, indent, name);
return;
}
// ***********************************************************************
//
// Works if the GeOb is a tuple of points from a cartesian product
// affine space (A x A x A x ...)
//
Simplex::Simplex(char* namein, GeOb &t)
{
static char* sig = "Simplex::Simplex(char*, GeOb&)";
int i;
int tuplesize;
GeOb new_obj;
SpaceType tupletype;
// Set the type:
type = SIMPLEX;
holds = SIMPLEX;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
// Check to see if the space is a CP Affine space:
tuplesize = (t.SpaceOf()).CPSpaceSize();
tupletype = (t.SpaceOf()).Holds();
if ((tupletype != AFF_SPACE) || (tuplesize == 1)) {
errh.ErrorExit(sig, "Space is not a cartesian product affine space", t);
}
s = t.SpaceOf()[0];
tostdb = Matrix(s.MatrixSize());
// Since the point tuple is from an affine space, each tuple element
// is an affine point. We just need to check that the tuple is the
// correct size and that all points are from the same space as we build
// the matrix:
if (tuplesize != s.MatrixSize()) {
errh.ErrorExit(sig, "Tuple is not the correct size for the space", t, s);
}
for (i = 0; i < s.MatrixSize(); i++) {
new_obj = t[i];
if (new_obj.SpaceOf() != s) {
errh.ErrorExit(sig, "Tuple element not from correct space", new_obj, s);
}
tostdb[i] = (new_obj.MatrixRep())[0];
}
// Make sure that the simplex is valid by checking that all the
// points are independent. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Tuple elements are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
//
// Builds a Simplex in the specified space:
//
Simplex::Simplex(char* namein, ASpace& ins, GeObList& t)
{
static char* sig = "Simplex::Simplex(char*, Space&, GeObList&)";
int i;
Boolean cannot_map = FALSE;
GeOb new_obj;
// Set the type:
type = SIMPLEX;
holds = SIMPLEX;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
s = ins;
tostdb = Matrix(s.MatrixSize());
// Check to see if we have the correct number of points, and build
// the basis matrix using the points:
if (t.Length() != s.MatrixSize()) {
errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
}
for (i = 0; i < s.MatrixSize(); i++) {
cannot_map = FALSE;
if (t[i].CanMapTo(AFF_POINT)) {
new_obj = t[i].MapTo(AFF_POINT);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else {
tostdb[i] = (new_obj.MatrixRep())[0];
}
} else {
cannot_map = TRUE;
}
if (cannot_map) {
errh.ErrorExit(sig, "Object cannot be mapped to specified space",
t[i], s);
}
}
// Make sure that the simplex is valid by checking that all the
// points are independent. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Points are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
//
// Used to cast a general basis down to a Simplex. This is also the
// way to change a Frame into a Simplex
//
Simplex::Simplex(Basis &f) : (f)
{
static char* sig = "Simplex::Simplex(Basis&)";
int i;
type = SIMPLEX;
// Conversion of Frames to Simplices:
if (holds == FRAME) {
holds = SIMPLEX;
// Add a prefix to the name:
char buf[MAX_NAMESIZE + 13];
char *smplxfrom = "Simplex from ";
strncpy(buf, smplxfrom, 13);
this->Name(buf + 13);
strncpy(name, buf, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
// Need to build the new matrix by taking the frame origin and
// adding it to all the vectors to make points. The origin stays
// as the last row of the target matrix.
for (i = 0; i < (s.MatrixSize() - 1); i++) {
tostdb[i] = (tostdb[s.MatrixSize() - 1] + tostdb[i])[0];
}
// Make sure that the simplex is valid. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Points are not independent", *this);
}
fromstdb = Inverse(tostdb);
} else if (holds != SIMPLEX) {
errh.ErrorExit(sig, "Attempted to cast a non-Simplex to a Simplex", f);
}
}
// ***********************************************************************
// ***********************************************************************
//
// Frame Class
//
// ***********************************************************************
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast a general basis down to a Frame.
//
Frame::Frame(Basis &f) : (f)
{
type = FRAME;
if ((holds != FRAME) && (holds != NULL_BASIS)) {
errh.ErrorExit("Frame::Frame(Basis &)",
"Attempted to cast a non-frame to a frame", f);
}
}
// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment:
//
Frame& Frame::operator=(Frame &f)
{
holds = f.holds;
s = f.s;
f.Name(name);
tostdb = f.tostdb;
fromstdb = f.fromstdb;
return (*this);
}
// ***********************************************************************
//
// Output for debugging:
//
void Frame::debug_out(ostream &c, int indent)
{
static char* name = "Frame object\n";
this->data_out(c, indent, name);
return;
}
// ***********************************************************************
//
// Builds a Frame in the specified space:
//
Frame::Frame(char* namein, ASpace& ins, GeObList& t)
{
static char* sig = "Frame::Frame(char*, Space&, GeObList&)";
int i;
Boolean cannot_map = FALSE;
GeOb new_obj;
// Set the type:
type = FRAME;
holds = FRAME;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
s = ins;
tostdb = Matrix(s.MatrixSize());
// Check to see if we have the correct number of objects, and build
// the basis matrix using the objects:
if (t.Length() != s.MatrixSize()) {
errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
}
// Map all but the last object to vectors:
for (i = 0; (i < (s.MatrixSize() - 1)) && (cannot_map == FALSE); i++) {
if (t[i].CanMapTo(AFF_VECTOR)) {
new_obj = t[i].MapTo(AFF_VECTOR);
if (new_obj.SpaceOf() != s.GetSpace(TANGENT)) {
cannot_map = TRUE;
} else {
tostdb[i] = (new_obj.MatrixRep() * s.HoistTangent())[0];
}
} else {
cannot_map = TRUE;
}
}
// Map last object to a point:
if (t[s.MatrixSize() - 1].CanMapTo(AFF_POINT)) {
new_obj = t[s.MatrixSize() - 1].MapTo(AFF_POINT);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else {
tostdb[s.MatrixSize() - 1] = (new_obj.MatrixRep())[0];
}
} else {
cannot_map = TRUE;
}
if (cannot_map) {
errh.ErrorExit(sig, "Object cannot be mapped to specified space", t, s);
}
// Make sure that the frame is valid by checking that all the
// vectors are independent. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Vectors are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
//
// Used to convert a Simplex into a Frame by selecting a
// point to serve as the origin:
//
Frame::Frame(Simplex& b, int n)
{
static char* sig = "Frame::Frame(Simplex&, int)";
int i;
int slot;
// Set the type, name, and space:
type = FRAME;
holds = FRAME;
s = b.SpaceOf();
tostdb = Matrix(s.MatrixSize());
// Prefix the name before copying:
char buf[MAX_NAMESIZE + 11];
char *framefrom = "Frame from ";
strncpy(buf, framefrom, 11);
b.Name(buf + 11);
strncpy(name, buf, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
// Check to see if the specified n is legal:
if ((n < 0) || (n >= s.MatrixSize())) {
errh.ErrorExit(sig, "Specified n is out of range",
ErrVal("n = ", n), b);
}
// Need to build the new matrix by taking the chosen point and
// subtracting it from all the the other points to make vectors.
// The chosen point is put in the last row of the target matrix.
// Map all but the last object to vectors:
for (slot = 0, i = 0; i < s.MatrixSize(); i++) {
if (i == n) {
tostdb[s.MatrixSize() - 1] = (b.Tostdb())[n];
} else {
tostdb[slot++] = ((b.Tostdb())[i] - (b.Tostdb())[n])[0];
}
}
// Make sure that the frame is valid. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Vectors are not independent", *this);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
// ***********************************************************************
//
// VBasis Class
//
// ***********************************************************************
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast a general basis down to a VBasis.
//
VBasis::VBasis(Basis& f) : (f)
{
type = VBASIS;
if ((holds != VBASIS) && (holds != NULL_BASIS)) {
errh.ErrorExit("VBasis::VBasis(Basis&)",
"Attempted to cast a non-VBasis to a VBasis", f);
}
}
// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment:
//
VBasis& VBasis::operator=(VBasis& b)
{
holds = b.holds;
s = b.s;
b.Name(name);
tostdb = b.tostdb;
fromstdb = b.fromstdb;
return (*this);
}
// ***********************************************************************
//
// Output for debugging:
//
void VBasis::debug_out(ostream &c, int indent)
{
static char* name = "VBasis Object\n";
this->data_out(c, indent, name);
return;
}
// ***********************************************************************
//
// Works if the GeOb is a tuple of vectors from a cartesian product
// vector space.
//
VBasis::VBasis(char* namein, GeOb& t)
{
static char* sig = "VBasis::VBasis(char*, GeOb&)";
int i;
int tuplesize;
GeOb new_obj;
SpaceType tupletype;
Space firstspace;
GeObType targettype;
Boolean cannot_map = FALSE;
// Set the type:
type = VBASIS;
holds = VBASIS;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
// Check to see if the space is a CP Vector space:
tuplesize = (t.SpaceOf()).CPSpaceSize();
tupletype = (t.SpaceOf()).Holds();
if ((tupletype != VEC_SPACE) || (tuplesize == 1)) {
errh.ErrorExit(sig, "Space is not a cartesian product vector space", t);
}
// Figure out the space to build the basis in. The vector tuple could
// be a mixture of vectors from a linearization/tangent space pair.
// Which space to use depends on the size of the tuple and the dimensions
// of the spaces:
firstspace = t.SpaceOf()[0];
if (firstspace.Dim() == tuplesize) {
s = firstspace;
} else if (firstspace.Dim() == tuplesize + 1) {
s = firstspace.GetSpace(TANGENT);
} else if (firstspace.Dim() == tuplesize - 1) {
s = firstspace.GetSpace(LINEARIZATION);
} else {
errh.ErrorExit(sig, "Tuple is not the correct size for the space",
t, firstspace);
}
targettype = s.NativeType();
tostdb = Matrix(s.MatrixSize());
// Now build up the matrix by mapping the vectors into the desired space:
for (i = 0; i < s.MatrixSize(); i++) {
cannot_map = FALSE;
if (t[i].CanMapTo(targettype)) {
new_obj = t[i].MapTo(targettype);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else {
tostdb[i] = (new_obj.MatrixRep())[0];
}
} else {
cannot_map = TRUE;
}
if (cannot_map) {
errh.ErrorExit(sig, "Tuple element cannot be mapped to specified space",
t[i], s);
}
}
// Make sure that the basis is valid by checking that all the
// vectors are independent. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Vectors are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
//
// Builds a VBasis in the specified space:
//
VBasis::VBasis(char* namein, VSpace& ins, GeObList& t)
{
static char* sig = "VBasis::VBasis(char*, Space&, GeObList&)";
int i;
Boolean cannot_map = FALSE;
GeOb new_obj;
GeObType targettype;
// Set the type:
type = VBASIS;
holds = VBASIS;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
s = ins;
tostdb = Matrix(s.MatrixSize());
// Check to see if we have the correct number of vectors, and build
// the basis matrix using the vectors:
if (t.Length() != s.MatrixSize()) {
errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
}
// Map list elements to correct types:
targettype = s.NativeType();
for (i = 0; i < s.MatrixSize(); i++) {
cannot_map = FALSE;
if (t[i].CanMapTo(targettype)) {
new_obj = t[i].MapTo(targettype);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else {
tostdb[i] = (new_obj.MatrixRep())[0];
}
} else {
cannot_map = TRUE;
}
if (cannot_map) {
errh.ErrorExit(sig, "Object cannot be mapped to specified space",
t[i], s);
}
}
// Make sure that the basis is valid by checking that all the
// vectors are independent. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Vectors are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************
// ***********************************************************************
//
// HFrame Class
//
// ***********************************************************************
// ***********************************************************************
//
// Public member functions
//
// ***********************************************************************
//
// Used to cast a general basis down to an HFrame.
//
HFrame::HFrame(Basis& f) : (f)
{
type = HFRAME;
if ((holds != HFRAME) && (holds != NULL_BASIS)) {
errh.ErrorExit("HFrame::HFrame(Basis &)",
"Attempted to cast a non-HFrame to an HFrame", f);
}
}
// ***********************************************************************
//
// Need to do memberwise assignment, since Matrix class members need
// to do memberwise assignment:
//
HFrame& HFrame::operator=(HFrame& b)
{
holds = b.holds;
s = b.s;
b.Name(name);
tostdb = b.tostdb;
fromstdb = b.fromstdb;
return (*this);
}
// ***********************************************************************
//
// Output for debugging:
//
void HFrame::debug_out(ostream &c, int indent)
{
static char* name = "HFrame Object\n";
this->data_out(c, indent, name);
return;
}
// ***********************************************************************
//
// Builds a Projective Frame in the specified space:
//
HFrame::HFrame(char* namein, PSpace& ins, GeObList& t)
{
static char* sig = "HFrame::HFrame(char*, Space&, GeObList&)";
int i;
Boolean cannot_map = FALSE;
GeOb new_obj;
Matrix unit_pt, tempinv, coeff;
// Set the type:
type = HFRAME;
holds = HFRAME;
// Local copy of the debug name:
strncpy(name, namein, MAX_NAMESIZE - 1);
name[MAX_NAMESIZE - 1] = '\0';
s = ins;
tostdb = Matrix(s.MatrixSize());
// Check to see if we have the correct number of points, and build
// up a temporary matrix and a unit point matrix using the points.
if (t.Length() != (s.MatrixSize() + 1)) {
errh.ErrorExit(sig, "Incorrect number of objects specified", t, ins);
}
Matrix tempmtx(s.MatrixSize());
for (i = 0; (i <= s.MatrixSize()) && (cannot_map == FALSE); i++) {
cannot_map = FALSE;
if (t[i].CanMapTo(PROJ_POINT)) {
new_obj = t[i].MapTo(PROJ_POINT);
if (new_obj.SpaceOf() != s) {
cannot_map = TRUE;
} else if (i == s.MatrixSize()) {
unit_pt = new_obj.MatrixRep();
} else {
tempmtx[i] = (new_obj.MatrixRep())[0];
}
} else {
cannot_map = TRUE;
}
}
if (cannot_map) {
errh.ErrorExit(sig, "Object cannot be mapped to specified space",
t[i], s);
}
// We need to scale the rows of the temporary matrix so that the
// unit point has homogenous coordinates (1, 1, 1, ...). Do this
// by solving a system of equations. Invert the temporary matrix:
if (fabs(Det(tempmtx)) < EPSILON) {
errh.ErrorExit(sig, "Points are not independent", t);
}
tempinv = Inverse(tempmtx);
// Solve for the row coefficients. If any coefficient = 0, the
// unit point was not in general position. Otherwise, multiply
// the matrix rows by the coefficients:
coeff = unit_pt * tempinv;
for (i = 0; i < s.MatrixSize(); i++) {
if (fabs(coeff[0][i]) < EPSILON) {
errh.ErrorExit(sig, "Unit point not in general position", t);
} else {
tostdb[i] = (coeff[0][i] * tempmtx[i])[0];
}
}
// Make sure that the hframe is valid by checking again that the
// matrix is invertible. If it is OK, invert it:
if (fabs(Det(tostdb)) < EPSILON) {
errh.ErrorExit(sig, "Points are not independent", t);
}
fromstdb = Inverse(tostdb);
}
// ***********************************************************************